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^ ' Abstract 



We have shown in previous work that statistical inference for cooper- 
ative sequential adsorption model can be based on maximum likelihood 
estimation. In this paper we continue this research and establish asymp- 
totic normality of the maximum likelihood estimator in thermodynamic 
limit. We also perform and discuss some numerical simulations of the 



ijD ' model. 

cn '. 

^T) . Keywords: cooperative sequential adsorption, time series of spa- 

CN I tial locations, spatial random growth, maximum likelihood estimation, 

lO ' asymptotic normality, Fisher information, martingale, thermodynamic 

Q ! limit 

2000 MSC: Primary 62M30, 60K35; Secondary 60D05 



X ■ 1 Introduction 



This paper continues the research started in [TO], where properties of maxi- 
mum hkehhood estimator for cooperative sequential adsorption model (CSA) 
were studied. CSA is a probabilistic model motivated by adsorption processes 
in physics and chemistry ([5]). The main peculiarity of adsorption processes 
is that adsorbed particles change adsorption properties of the material. For 
instance, the subsequent particles might be more likely to be adsorbed around 
locations of previously adsorbed particles. In other words, the adsorption pro- 
cess might accelerate as the surface gets saturated. In the opposite scenario 
adsorbed particles inhibit adsorption of subsequent particles, so that the ad- 
sorption process slows down. 
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Mathematically CSA is formulated as a random sequential allocation of 
points in a bounded region of space (the observation window). The result of 
CSA dynamics is a sequential point pattern which seems to be of great interest 
in many apphcations. It should be noted that CSA can produce a large variety 
of aggregated point patterns (see, e.g., the images throughout the paper). 

It was first noticed by physicists (e.g., see [S], p. 1285) that this type of 
model can be used for modelling the spatial-temporal processes similar to the 
irreversible spread of disease or epidemics. This idea is developed further in 
[To] where use of CSA for modelling time series of spatial locations is discussed. 

Biological growth was mentioned in [5] as another potential application of 
the adsorption models. These ideas have been recently supported by both ex- 
perimental and simulation studies of keratin filament (KF) network formation 
in biology. KF networks are part of the cell cytoskeleton and they determine 
the shape and biophysical properties of the cells. Loosely speaking, the KF is 
an aggregated spatial structure formed by a union of curved finite segments 
(fibres). Experimental results ([Hj) and simulation studies ([!]) suggest that 
the KF can be thought as a result of a sequential spatial growth process with 
self-organising properties. It is also argued in [^j (see also references therein), 
that self-organizing processes combined with simple physical constraints seem 
to have key roles in controlling organelle size, number, shape and position, 
and these factors then combine to produce the overall cell architecture. CSA 
seems to be useful for modelling spatial random growth with self-organising 
properties. 

The variant of CSA under consideration here is easy to parameterise. Sta- 
tistical inference for the model parameters developed in [10] was based on 
maximum likelihood estimation (MLE). It was shown in [10] that maximum 
likelihood estimator exists uniquely. Moreover, it was proved that the maxi- 
mum likelihood estimator is consistent in the thermodynamic limit. The ther- 
modynamic limit means that the observation window expands to the whole 
space and the number of allocated points grows linearly in the volume of the 
window. The main result of the present paper is asymptotic normality of 
maximum likelihood estimator in the same limit. 

The study of statistical properties of MLE in both [10] and this paper is 
essentially based on the fact that the model likelihood depends on the point 
configuration via statistics with a certain special structure, allowing us to 
apply the limit theory for random sequential packing and deposition (see, e.g., 

2 CSA as a generalisation of random sequen- 
tial adsorption 

The adsorption model most commonly studied in the physics literature is ran- 
dom sequential adsorption (RSA). Mathematically RSA is formulated as the 



following packing model. Consider a bounded region D of Euclidean space 
R (modelling the adsorbing material) and a sequence of independent points 
Yi, ^2, • • • , (modelling the particles) sequentially arriving in D at random. An 
arriving point is accepted with probability 1, if the ball of a certain fixed radius 
R (interaction radius) centered at the point does not cover any of previously 
accepted (adsorbed) points; otherwise the point is rejected. 

RSA with interaction radius R is nothing else but the ci- dimensional ver- 
sion of the classical car parking model [11], where a "car" is modelled by a 
ball of radius R/2. Clearly the distance between any two points in a RSA 
point pattern cannot be less than the interaction radius R. Therefore RSA 
generates only regular point patterns which are similar to the right one in Fig- 
ure [2], and never generates point patterns similar to the left one in Figure [2J 
However, RSA can be easily generalised in order to generate aggregated point 
patterns. To do so, we allow neighbours. That is, we let an arriving point 
be accepted with a certain conditional probability, even if a ball of radius R 
centred at the point covers some of the previously accepted points. In general, 
the acceptance probability can depend on the spatial configuration formed 
by previously accepted points. We study the model in which the acceptance 
probability depends on the number of neighbours. 

More precisely, fix a sequence of non- negative numbers /3o, /3i, ..., such that 
/3o > 0. Given a sequence of accepted points X{k) = {Xi, . . . ,Xk) (-^(0) = 
0), let the next uniform arrival Y be accepted with conditional probability 
proportional to /3j, if the number of neighbours of Y among Xi, . . . ,Xk is 
equal to z > 0. If /3q > and /3fc = 0, A; > 1, then this model is RSA. 

This CSA model can be regarded as a continuous version of the lattice 
model (i.e. where Z^ is a subset of lattice Z ) known as monomer filling 
with nearest-neighbour cooperative effects. CSA in this particular form was 
formulated for the first time in [12] , where its asymptotic study was undertaken 
under certain assumptions. In what follows we denote by CSA the adsorption 
model of this type. 

CSA can be used for modelling both clustered and regular point patterns. 
A large variety of aggregated point patterns can be generated by modulating 
the model parameters. For instance, the left image in Figure [21 containing 
1000 points, is generated by CSA with parameters R = 0.01, /3o = l,/3i = 
1000, /32 = 10000, /3fc = 0, A; > 3. The right image (containing 500 points) 
is a typical regular image produced by RSA (here the interaction radius is 
i? = 0.03). 

3 Notation and assumptions 

Let D he a, convex compact subset of R , i? be a positive constant, and 
{/3fc, k > 0} he a sequence of non- negative numbers. For any point x G R 
and any finite sequence y = (i/i, . . . , y„), n > 1, of points in R*^, we denote 
by z^(x,y) the number of points yi in the sequence y, such that the distance 
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between x and |/j is not greater than R. By definition u{x, 0) = 0. 

Let X{i) = {Xi, . . . , Xi)^ Xi E IV^, z = 1, . . . , £ be a vector of first i random 
points sequentially generated by CSA. CSA dynamics goes as follows. Given 
a sequence of accepted points X{k) = {Xi, . . . ,Xk) (which can be empty, i.e. 
A; = 0) a new point Y, uniformly distributed in D, is accepted with probability 
proportional to /3u{Y,x{k)) and rejected otherwise. If Y is accepted, then we set 
Xk^i = Y and X{k + 1) = (Xi, . . . ,Xfc,Xfc+i). The conditional probability 
density function of the next accepted point X^+i is 



i'k+iix) 



/3. 



(x,X(k)) 



Jijl3uiy,xik))dy 



X e D. 



It is easy to see that the sequence of accepted points is an embedded Markov 
chain for a continuous time spatial birth process x(t) E D, t > 0, specified by 
the following birth rates. If the process state at time t > is x, then the birth 
rate at point x G -D is /3u(x,x.), the total birth rate is 



ax 



Pv(x,-x.)0''X-i 



D 



and the waiting time until the next process jump is an exponential random 
variable with mean a~^(x). 

As in [To], we assume throughout that 

• there is a finite number of positive /3's, that is /3o > 0, . . . , (ij^ > and 
/3fc = 0, for k > N + 1, for some N > 1, where the number N can be 
unknown. 



/So 



• the interaction radius i? is a fixed and known constant. 

It is easy to see that the joint probabihty density Y[k=i'^k{xk) of the first 
accepted points can be written as follows: 






Pe,i3,Dixi, ...,xe)- ^^ ^^l "^ r l(7Vfem)<m' (2) 

where 



N{x{t))= max u{xi,x{i — 1)), (3) 

and 

i 
tkW)) = ^ l{^(^,,3;(i_i))=fc}, k = 0,...,N. (4) 

j=i 

where we denoted for short x{k) = {xi, . . . , Xfc), fc > 1, and x(0) = for k = 0. 
Remark. It should be noticed that we do not completely recover the pa- 
rameters of the spatial birth process. In the present setting we do statistical 
inference only for the embedded Markov chain, which distribution is completely 
specified by the ratios Pi/Po^ 1, . . . ,N and the interaction radius. As a result, 
one can forecast the probability distribution of the next accepted point, but not 
the waiting time until the next acceptance event. 

As in [TU], let Di be the unit cube centred at the origin and consider a 
sequence of rescaled domains 

D^ = m^/'^Di, m G Z+. 

Fix {im, "m > 1} an arbitrary monotonically increasing sequence of positive 
numbers, where im stands for the number of observed points in the domain 

Dm. 

Assumption 1 The number of observed points is asymptotically linear in m, 
that is 

lim f^) =yuG(0,M, 

where O^d is the jamming density (fl^). 
Define 

S.m := {x{im) = {xi, . . .,xej, Xi e Dm : NWm)) < N}. 

Given parameters /3 = (/3i, . . . , Pn) consider a probability measure Pm,/3 on Sm 
specified by the probability density (|2]) with i = im and D = Dm- Expectation 
with respect to this measure is denoted by Em,/3- We assume that /3 E B, 
where B is an open subset of R^, such that B C R^. The true parameter 

is denoted by /3^°^ = f /3| , . • • ^ Pn )■ Also, we denote for short Pm = Pm,/3(o) 

and E^^ = E^_^(o) . 



4 The results 

Given m assume im^ '^ and consider the log likelihood function 

Lm{X-{iJ, (5) = \og{pe^,^,n^{Xr, ..., X™ )), (5) 

where X"^{im) is the vector of observed points in D^n- Given observation 
X"^{im) we define the maximum likelihood estimators 

of parameters /3(°) = (/3| % . . . , /3]y ) as maximizers of function Lm{X^{lrn)i P) 
and which can be found as a solution of the following system of MLE equations 

^ r,n ' =0. J = l,...,iV. 6 

The following two statements were proved in ^U\ (see Theorem 2.2 and 
Lemma 5.2, part 2), respectively in [TO]). 

Lemma 4.1 Under Assumption{l\ with Pm— probability tending to I as m ^ 
oo there exists a unique positive solution (/3i,m, • • • , /Sat,™.) of the likelihood equa- 
tions and 

(/3l,m, • • • , (3N,m) — ^ (/3l , • • • , /^AT ) 

in Pm — probability as m -^ oo. 
Lemma 4.2 Consider the matrix 

MX (U,/3).--^ g^^^ )^^^^^. 

There is a family of N x N real matrices J{(3,fj,), defined for (3 E B and 
fi G (0,^oo); such that under Assumption^ 



m 

in Pm — probability as m ^ oo for any (3 E B. Moreover, the limit matrix 
evaluated at (3 = (3^°\ i.e. 

jW(/i) = j(/3W,^) (7) 

is positive definite. Finally, if (3{m) is a random B-valued sequence converging 
in probability to (3^^^^ as m —> oo, then 

m 
in Pm — probability as m ^ oo. 



The last part of Lemma [4.21 is not included in Lemma 5.2(2) of [TU], but can 
be proved in the same manner as that result. 

In Section[7]we give extended study of the structure of the limit information 
matrix. 



Theorem 4.1 Under Assumption^ the model score function 



(0), ldLm{X^^{im),P) aL„(X-(C),/3) 



13=13(0) 



converges in distribution as m —)■ oo to a Gaussian vector with mean zero and 
covariance matrix J*-°^(/i). 

Theorem 14.11 is proved in Section 16.21 The following theorem states that the 
MLE is asymptotically normal. This is the main result of the paper. 

Theorem 4.2 Under Assumption^ 

in distribution as m ^ oo, where Af iO, {J^^^) (/i) ) is the Gaussian vector 
with zero mean and with the covariance matrix (j(°)(/i)) 

Theorem 14.21 provides asymptotic justification for creating confidence intervals 
based on the normal distribution, as we do in the example in Section [3 

5 The model likelihood 

In this section we introduce more notation and recall some other facts from 
[To] which will be used in Section |6l 

Let X^{im) = (X^, . . . ,X^) be the sequence of observed points X™ in 
Dm- Denote 

t]^, = tj{X^ik)) 0<k<im-l,j = l,...,N, (9) 

where tj, j = 1, . . . ,N, are statistics defined by equation (jl]), and denote 

^Tk = rrfc(^"(^)) = / Uu:uiu,X^{k))=j}du, 0<k<im-l,J>0, (10) 

Dm 

note that T]^,^{X'^{k)) = ioi k < j and that T}^Q{X"'{k)) is equal to m. 

In terms oft— and F— statistics, using ([2]), (jl]) and flTOj) the model likelihood 
can be rewritten as follows 

L^(x'"(c), P) = iog(p,„,^,z,^(xr, . . . , X™ )) 

N 



JV tm / p 

J]tfc(X™(C)) log(/3fc) - J^log ( / f3uiu,x^ik-i))du 

fe=l fe=l y-'Dm 

N im / N \ 

E C^™ iog(/5'^) - E log r™,_, + E /3.r,>-i • (ii: 

fc=l k=l \ j=l J 



Thus the log hkehhood function depends on the observed point configura- 
tion only through t— statistics tk and F— statistics V^,^. 
Theorem 2.2 in [10] says that if 

lim(£^/m)=^G(0,^oo(/3)), 

then as m — )■ oo we have for any (5 ^ B that 



m 
and 



^■'^-'"4p,(/i,/3),j = l,...,iV, (12) 



^■'^-''"4 7,(/^,/3),J = 0,...,iV, (13) 



m 



where the functions {pj (/i, /3) ,/i G (0, 6'oo(/3)), 1 < j < A^ and (7j (//, /3) ,/U G 
(0, 6'oo(/3)), < j < A^ are strictly positive and continuous in /x, and are related 
by the following integral equation 

J 7o(A,/3) + )_^.^^/3i7i(A,/3) 
for any < /i < 0oo{f3). 

6 Proofs 

6.1 Proof of Theorem 11:21 



Given Theorem I4.H the proof of Theorem \A.2\ although new to this particular 
model, runs along standard lines (see e.g. [1], or Theorem 1 of [2]), and we 
give just a sketch. 

Choose S > such that the ball of radius S centered at /3^'^^ is contained in 
B. By consistency of the maximum likelihood estimator P{X"^{im)) (Lemma 
14. ip . we have that 

|^(X™(U)-/3W|<5, 

with probability Pm close to 1 if m is large enough. With di denoting differ- 
entiation with respect to the ith component of /3, we make a Taylor expansion 
of 9,(L^(X™(£^),/3) about /3W: 

= d,LUX"\im),Ax"'{iJ)) = 9,L„(X"^(£„),/3(°)) 

N 

+ ^9J(X'"(C) J)(^(x™(0) - /3^°^)i, 



where (3 lies on the hne segment from (3^'^'^ to /3(X™(£m)). Rewriting this 
expression, we obtain 

In the left hand expression /3 depends on i but converges in probability to 
/S^"-* as n — )■ oo by Lemma [4. II By Lemma \A.2\ for each (i, j) the first factor 
inside the sum converges in probability to Jl: {jj). Observing that Theorem 
14.11 applies to the right hand side, we can complete the proof by applying 
Lemma 6.4.1 of [1]. 

6.2 Proof of Theorem liJi 

Let J-"]™ = cr{Xj^, . . . , X™} be the a— algebra generated by the first j points 
observed in Dm- Asymptotic normality of the score function is essentially 
based on the following fact. Namely, for any fc = 1, . . . , A^, the triangle array 



tm,/3 



dPk 



-j,^ra)\j,im)\^ , m > 2, (15) 



J=l 



is a zero-mean square integrable martingale array. Indeed, by the representa- 
tion (ITTD. 



for j = 1, . . . , A^. Introducing the following quantities 

^k',i = l{i^(X™,X'"(i-l))=A:}, k = 0,. .. ,N, i = 1,. . .,irn, (17) 

allows to rewrite equation (jl]) as follows: 

tfc(X'"(£„)) = ^eM, A; = 0,...,Ar. 

Denote for short 

d^T- p / ^^"^ I Tr("^) \ 

It is easy to see that 

4- = -^^^ , k = l,...,N, t = l,...,im- (18) 

By using notation 

Cfc,j = 'W' \^k,i ~ ik,i) (19) 

Pk 



equation (1161) can now be rewritten as follows: 

dL^X-iiJ, /^) ^ 1 ^ ^^rn^ _ -^rn^ ^ ^ ^^ fc = 1, . . . , iV. (20) 

Therefore the triangle array flTSl) is a zero-mean square integrable martingale 
array with differences given by equation flT9|) . This implies that for any real 
vector a = (ai, . . . , a^Y 



l_^ gL^(X"-(C),/3) 



k=l 






is a zero-mean square integrable martingale array. 

By the Cramer- Wold device (see for example [3]), Theorem 14. 1 1 follows from 
the following fact. 

Lemma 6.1 Under AssumptionUl for any real vector a = (ai, . . . , a^)'^ , 
1 ^ dLUX-{U,/3('^) ^ 

m distribution as m ^ oo, where 



-' "^J(°V)J 



cr., = a J^ M u a, 



J'^°^(/i) zs i/ie matrix defined by equation ([^ and A/'(0,cra) is ^/ie Gaussian 
vector with zero mean and variance a^. 

In proving Lemma 16.11 we shall repeatedly use the following fact which is 
simple enough for us to omit its proof. 

Proposition 6.1 Let C,n, n > 1, and r]n, n > 1, be two sequences of random 
variables, C > 0,a and b be some constants. Suppose that |^„| < C, \rin\ < C, 
^n -^ a in probability as n —>■ oo and E(?7„) —)■ 6 as n —)■ oo . Then E(^„?7„) — )■ ab 
as ra —)■ oo. 

Proof of Lemma \6.1\ By (12U]) , for any 13 & B we have 

^A gL^(X-(C),/3) _ 1 y^ m 

where 

N 



Y.^kcz (21) 



fc=i 



and Cfc^j are the quantities defined by equation (IT^ . It is easy to see that 

— ^= max InJ"! < --= max ( tt I — ^ 0, as m — )■ oo, (22) 



and 

E^.5 ( max (ryrf)< — 
m 

By Propositions 16.21 and 16.31 below, we also have under Assumption [1] that 



— Em/3 (max(?7,'") ) < max ( -— | — t- 0, as m — > oo. (23) 

m '^ \ i ^ " ' J m k=i,...,N \f3kj 



^Em^^a^j(°)(/.)a (24) 



m 

i=2 



in Pm— probability as ?Ti — ?■ oo. 

Using fl2^ . fl2^ and fl2^ . we can then apply the central limit theorem for 
martingale difference arrays (Theorem (2.3) of [7]) to complete the proof of 
Lemma 16.11 



Proposition 6.2 Under Assumption\^ 

m— i>oo rn 



li--EEL°U(Cn=a^^^°H/^)a. (25) 

n— i-no m ' ' ^ ' 



i=2 



Proof of Proposition \6.2[ It was shown in Section 6.2 of [10] that the limit 
of the scaled Hessian in Lemma 14.21 evaluated at the true parameter has the 
following integral representation 

j(°)(/x) = j(/3(°),/x)= fQ^'^\X)dX, (26) 



where 

where 5ij is the Kroneker symbol, 7] (A) = 7j(A, /S^^-*), j = 0, . . . , A^ (7— functions 
are defined by ( IT3l) ) and 

AT 

Z(/3,A)=7r(A) + E/5^^1°^(^)- (28) 

Let us show that ii i = im is such that i/m — )■ A G (0, /i), as ?7i — )> 00, then 

Ei?((C)^)^a^gW(A)a 



as m — !• oo. Indeed, 

N 



^^ iivrr) = E «^«.EL°nQcrj 



k,j=i 

N 

/q(0)/q{0) "^ ^^'^'='* ^k,i)\.^j,i ^j,i)l 



N 

Qkttj 



k,j- 



fi ^r^r 



E ^ J p(0) ( C'nT' Ci^ C"^ C iT^ P"^P "^ -\- ^ "^/^ "^\ 
0(0)0(0) 'Ti \^k,i^j,i ^k,i>:j,i Sj,jSfc,i "T Sfc,jSj,iJ 

k,j=l Pk Pj 



Notice that 

Em {<.k,i^j,i) = E^ {^k,i)^kj = E^ [ik,i)^kj-, 

where (5jj is the Kroneker symbol and 

m \^k,i^j,i J m \S)j,iSfc,i/ ra \^k,i>'j,i J ' 

Therefore 

N 
Em ((^I") ) = 2^ (O)o(O) Fm {ik^d^kj - E^ (Cfc™^i^)J 

By ([18]) and ^, 



" r- _, + Ef., /3f r-_, ^ 7r (A) + E,t, /^iSf (A) ^ (/^^^^ A) 

(29) 
in Pm probabihty as i/m — )■ A, for any r = 0, . . . , A^. This fact along with 
Proposition 16.11 yield that 



) 0(0) L m \'^k,i)"kj i-m VS>fc,jSj,i 



p'^'Pf 



-)■ 



/3;.'"Z(/3(»),A) " Z'(;3(»).A) '^'■'^^' 



as z/m — )■ A. We can then complete the proof of Proposition 16.21 by applying 
the dominated convergence theorem to show the sum converges to the integral 
(see Section 5.2 of [lO] for a similar argument.) 



Proposition 6.3 Under AssumptionUl 

0, (30) 

where the expectation is taken with respect to measure Pm ■ 



hm i-Var( V(r/nM 



Proof. To simplify notation we assume in the proof that A^ = 1; modifications 
for the multivariate case are obvious. Also, for simplicity of notation, we 
omit the upper index in notation for r], ( and ^ variables. So, in the rest of 
the proof we denote /3 = (3i, a = a e R, rji = rj^, Q = (^}, ^i = ^^i, ^i = 

^]^, J-'j = J-"! . Besides, we write E instead of Em . 
It suffices to show that under Assumption [T] 

Cov(r/2,r/|)^0 (31) 

for any pair of sequences i = im and j = jm such that i j^ j and i/m — )• 
A', j/m — )■ A" as m — !■ oo, where A' can coincide with A". This suffices because 
the contribution from terms with i = j, divided by m^, is asymptotically 
negligible since the rji are uniformly bounded. 

Recall that rji = a(^j — ^i)//3, where ^j = E(^j|J^j„i). Therefore, we need 
to prove that 

Cov((e.-6)',(0-^.)')^0 (32) 

under the same assumptions about the index sequences. Assuming for definit- 
ness that i < j, we have the following identities 

C' = e^, (33) 

E(e.) = E(E(e.|J-.-i)) = E(e,), (34) 

E(a.) = E (6E (e.|-F,_i)) = E(e7), (35) 

E(/(e.,6,o)e,) = E (/(e.,e.,o)E(e,i-^,-i)) = E(/(e.,e.,o)o), m 

where f{^i,^i,C,j) is a polynomial function, e.g., ^i^f etc. (note that (15^ fails 
for i = j.) We can write Cov ((,^j — ^i)^, (^^ — ^jY) as a linear combination of 
terms of the form 

E(eTe- -^aer"^) - E(eTe- -^) H^U]-') (37) 

where p G {0, 1,2} and q G {0, 1,2}. As mentioned before (see display (!29|) ). 
we have that 

as r/m -^ A, and also, E{C,r) — > b{X) as n — > oo. Since C,i and ^^ are bounded, 
we have E{^f) -^ b^X') and (using ([35])) E(^ifi) -^ b\X'), while (using (03])) 
E(^f) — 7- &(A'), and likewise for j. Therefore 

E(efe-"'')E(eKr') ^ 6i+--(P'i)(A')6'+°"°(''''nA")- (39) 

But using (!36|) . ( 138|) and Proposition 16. II we also have 

E(a--%0) = E(a--"^|) ^ 6^+-''^(^'^)(A')6^(A"), (40) 

and using fl33|) for j, fl36|) . f l38|) and Proposition 16. II we also have 

E(eTe'"%') = E(eTe-"%) ^ fe^+'^^-^^^'^HAOKA")- (41) 

Combining ( 140|) and ( 14T|) shows that E(^f^^~^^J^^~'^) converges to the same 
limit as the expression in (15^ . Hence, each expression of the form in (13 7p 
tends to zero, and we have established ( 15^ . Hence, Proposition 16.31 is proved. 



7 Numerical example 

In this section we give a numerical example demonstrating that MLE is effec- 
tive in distinguishing between CSA's which might generate similar patterns. 

In [To] we briefly discussed difference between clustering effects produced by 
CSA determined by a set of increasing parameters /3 (the so-called Aarhenius 
rates, [5]), and determined by a set of flat rates (the so-called Eden rates, [5]). 
As before, we consider two single realizations of CSA. Six successive images for 
each of realization shown in Figures [T]l6l The interaction radius is i? = 0.02 in 
both cases. The left images have been generated by CSA with /3-parameters 
(3o = l,/3i = 300, /32 = 500, /3fc = 0, k > 3. The right images have been 
generated by CSA with /3-parameters f3o = 1, (3i = (32 = 100, f3k = Ok > 3. 
The first five pairs of images with first i = 200, 500, 1000, 2000 and i = 3000 
points respectively are shown in Figures [l]l5l The last pair of images shows 
the realisations at jamming, i.e., when there is no space left to accommodate a 
point. The left image contains i = 4407 and the right image contains i = 4416 
points. Can one tell apart these two sets of parameters given the series of 
images provided? 

The images with 200 points look similar and it seems plausible that they 
have been generated by the same CSA. In both cases new points tend to 
appear in the vicinity of existing points because of the choice of the parameters. 
Though clusters formed by a single point are noticeable on the right image and 
clusters seem to be more dense on the left one. 

The pair of subsequent images containing 500 points is shown in Figure [2l 
It is noticeable at both images there are almost no new clusters; the existing 
clusters keep growing and eventually start coalescing. Besides, it is slightly 
visible that the right pattern is more dispersed than the left one. All these 
effects are becoming more visible for the pair of images showing further evolu- 
tion and containing 1000 and 2000 points. These images are shown in Figures 
SandU 

The effects that have been just described are rather straightforward ana- 
logues of the phenomenon of "competition between the birth, growth and 
coalescence" ([3], p. 1307), which is well known for lattice CSA models. 

Though the main basic feature of both series of images, namely, clustering, 
is common to both choices of the parameters, the clustering effect is more 
visible in the images produced by the model with an increasing set of non-zero 
parameters (the sequence of left images). The clusters are more saturated in 
the left images, i.e. clustering is stronger. It seems that the right realisation 
spreads faster in comparison to the left one. This is called mild clustering; 
the distribution of points inside a cluster is more or less regular, since a new 
point distribution is uniform conditioned on being adsorbed in the vicinity of 
existing points. 

The difference between the strong and the mild clustering (corresponding 
to increasing and fiat sets of non-zero parameters respectively) observed in 
Figures [Jllll vanishes at the later stages of evolution, when it approaches jam- 



ming. It is quite difficult to distinguish by visual inspection the two sets of 
parameters given the pair of images shown in Figure O Note that both of these 
images are close to the corresponding jamming images shown in Figure [6l One 
might argue that these two realisations have been produced by the same model 
and the differences between them (observed at some intermideate images) can 
be attributed to variability of the samples. Numerical results given in Tables 
[T] and [2] show that MLE is an effective tool for parameter estimation. The 
tables contain MLE's for both sets of parameters along with corresponding 
approximate confidence bounds (any computed value is rounded to its near- 
est integer). The 95% confidence bounds are computed by formally assuming 
normality of /3. The variances of the estimates are approximated, as usual, by 
the corresponding diagonal elements of the matrix inverse to the observed in- 
formation matrix. The latter turned out to be non-degenerate for all observed 
images. The variances of the estimates decrease as the number of observed 
points increase. As a result, the confidence intervals become narrower. The 
tendency breaks down only for the rightmost entry of the bottom line in Table 
[U Perhaps this can be explained by the lack of accuracy of the computations 
(see the discussion of computational issues in [10]). The observed reduction 
of variances is intuitively expected, although the normality assumption in the 
unit volume cannot be based on our asymptotic results. This is in contrast to 
the limiting situation where the effect is clearly implied by the integral rep- 
resentation (126|) for the information matrix. The representation implies that 
the variance of the estimate Pi converges, as m — )■ oo and im/i^ ~^ f^, to 



J^gii{X)dX' 

where gu^X) > is the i—th eigenvalue of matrix Q^'^\X) in the representaion 
( 126|) . The preceeding display justifies "reduction of variances" effect, if the 
density of points, i.e. /i, increases. The lower bound for the variance of the 
estimate Pi is given by the same formula with yU = 6*00, where 6^0 is the jamming 
density ffTO]). 

Under certain assumptions normality of /9 in a fixed finite volume can pos- 
sibly be advocated as follows. Consider, for definiteness, the model in the unit 
volume and let the interaction radius be sufficiently small. This is the case in 
the simulated examples. If the interaction radius is sufficiently small, then the 
jamming density is high. In other words, a sufficiently large number of points 
can be accommodated. It was shown in Section I^T^ that the score function is a 
martingale sum containing i terms, where i is the number of observed points. 
Therefore, one might expect that if £ is sufficiently large (e.g., thousands), 
then the normal approximation starts working. 

Finally, it should be noted that MLEs effectively capture the correct mag- 
nitude of the parameters and this is why two considered sets of parameters 
in the example (producing sometimes quite similar images) can be effectively 
distinguished. For the sake of completeness, consider also the left image in 



Figured It has been generated by CSA with the interaction radius 0.01 and 
/3 -parameters /3o = 1.0, A = 1000.0, /^a = 10000.0, /3fc = 0.0, A; > 3. The im- 
age contain I = 1000 points, t— statistics are to = 23, ti = 149,^2 = 828. The 
MLE estimates for /3i and P2 are 1105.0 and 10510.0 respectively. 
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Figure I: £ = 200. Left: increasing rates, (^0,^1,^2) = (16,93,91). Right: fiat rates, 
(to,ti,t2) = (43,100,57). 
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Figure 2: i = 500. Left: increasing rates, (^0,^1,^2) = (25,233,242). Right: fiat rates, 
(to,ii,i2) = (62,272,f66). 




Figure 3: i = 1000. Left: increasing rates, (^0,^1,^2) = (34,434,532). Right: flat rates, 
(io,ii,i2) = (84,552,364). 
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Figure 4: £ = 2000. Left: increasing rates, (to,ii,i2) = (43,825,1132). Right: flat rates, 
(to,ii,i2) = (95,1048,857). 





Figure 5: ^ = 3000 Left: increasing rates, (^0,^1,^2) = (47,1190,1763). Right: flat rates, 
(io,ti,t2) = (106,1473,1421). 





Figure 6: Left: increasing rates,^ = 4407, {t^MM) = (48,1426,2933). Right: flat rates, 
i = 4416, (<o, ^1,^2) = (108, 1688, 2620). 



Table 1: MLE's for the left images in Figures [TM] 



^ = 200 


£ = 500 


i = 1000 


i = 2000 


i = 3000 


i = 4407 


/3i = 401 
(176, 626) 


(3^ = 377 
(214, 540) 


/3i = 334 
(213, 455) 


/3i = 320 
(218, 422) 


/3i = 318 
(223, 413) 


/3i = 323 
(226, 420) 


/32 = 695 
(298, 1091) 


[32 = 594 
(335, 853) 


[32 = 566 
(360, 772) 


[32 = 546 
(371, 721) 


[32 = 521 
(364, 678) 


[32 = 533 
(373, 693) 



Table 2: MLE's for the right images in Figm'es [T][6] 



£ = 200 


£ = 500 


£ = 1000 


£ = 2000 


£ = 3000 


£ = 4416 


/3i=89 
(55, 123) 


/3i = 98 
(69, 127) 


/3i=96 
(73, 119) 


/3i = 104 
(81, 127) 


A = 101 
(80, 122) 


/3i = 98 
(78, 118) 


/3i = 106 
(61, 151) 


/3i = 97 
(67, 127) 


/?2 = 88 
(66, 110) 


[32 = 100 
(78, 122) 


/32=99 
(78, 120) 


/32 = 98 
(78, 118) 



Appendix. On positive definiteness of the limit 
information matrix 

It is easy to see from equation (!26|) that positive definiteness of matrix Q*-°^(A) = 
Q (/3^°'', a) for any fixed A G (0, ^oo) implies positive definiteness of the limit 
matrix j(°^(/i). Positive definiteness of matrix Q^^\\) was shown in Lemma 
5.2 in [To]. Here we give another proof by studying the matrix structure in 
more detail. 

It can be seen from equation (l27j) that the matrix principal minor formed 
by the intersection of the first k rows and k columns is 



/^iv,fc(/3(°\A) 



/ 7l°>(A)(z(/3(°),A)-7r>(A)/3f)) 
/3f)z2(/3(0),A) 



7f'(A)7r(A) \ 

Z2(/3(0)^A) 



v 



7r(A)7f'(A) 



7l°>(A)(z(/3(0),A)-7l°'(A)/3r) 



Z2(^(o),A) 

It is easy to see that determinant of D^^^if^^^K A) is 

{ -\\k ^ 



/3^'"Z2(/3(0),A) 



/ 



U^^,fe(/3^"SA) 



Z2/=(/3(o),A)il ^, 



n^i^^-^(/5^°^^)^^i 



(0) 



where \Ak — Z((3^^\ \)Ek\ is determinant of matrix A^ — Z[(3^^\ X)Ek, where, 
in turn, matrix Ak is defined as follows 



A 



.= (/3r,...,/3r)(7r(A),...,7r(A))^ 



(42) 



and Ek is the k x k unit matrix. By definition, |/lfc — Z{(3^^\ \)Ek\ is the 
characteristic polynomial of A^ evaluated at point Z{(3^^\ A). It can be shown 
(we omit the proof) that if a, fe G C" are non-zero complex vectors, such that 
a^b 7^ 0, then a quadratic matrix M = ab^ has the only non-zero eigenvalue 
a^b of multiplicity 1, is the other matrix eigenvalue of multiplicity n — 1 and 
the matrix characteristic polynomial is 

\M-uE\ = {-iy'u''-^{u-a^b), ugC". 

Hence, 

(N ^ 

7r(A)+ E/5i^^^°'w 
i=k+l 

and 

The right side of the preceding display is positive because the functions 7,, i = 
1, . . . ,N are positive. Thus any principal minor of matrix (127|) is positive and 
by Sylvester criterion this matrix is positive definite. 
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